#ifndef __OPERATOR_H__
#define __OPERATOR_H__

#include "Poisson.h"
#include <stdexcept>

#define TEMPLATE template <typename Element>

template <typename Element = P1Element>
class RestrictionOperator
{
protected:
	size_t n;

public:
	RestrictionOperator(size_t _n);
	virtual ~RestrictionOperator() = default;
	virtual const Vector operator*(const Vector& x) = 0;
	size_t GetN() const;
};

TEMPLATE
RestrictionOperator<Element>::RestrictionOperator(size_t _n) : n(_n) {}

TEMPLATE
size_t RestrictionOperator<Element>::GetN() const
{
	return n;
}

template <typename Element = P1Element>
class FullWeighting : public RestrictionOperator<Element>
{
public:
	FullWeighting(size_t _n);
	~FullWeighting() = default;
	virtual const Vector operator*(const Vector& x);
};

TEMPLATE
FullWeighting<Element>::FullWeighting(size_t _n) : RestrictionOperator<Element>(_n) {}

TEMPLATE
const Vector FullWeighting<Element>::operator*(const Vector& x)
{
	size_t n = RestrictionOperator<Element>::n;
	Element element;
	if (x.rows() != element.GetNodeNum(n))
		throw std::length_error("FullWeighting: sizeof(x) error!");
	Vector y(element.GetNodeNum(n / 2));
	y.setZero(element.GetNodeNum(n / 2));
	size_t k0, k1, k2, k3, k4, k5, k6, k7, k8;
	double num = sqrt(element.GetNodeNum(n)) + 1.0;
	num = n / num;
	for (double i = 2 * num; i <= n - 2 * num; i += 2 * num) {
		for (double j = 2 * num; j <= n - 2 * num; j += 2 * num) {
			k0 = element.SerialNumOfIndex(i - num, j - num, n);	
			k1 = element.SerialNumOfIndex(i, j - num, n);
			k2 = element.SerialNumOfIndex(i + num, j - num, n);
			k3 = element.SerialNumOfIndex(i - num, j, n);
			k4 = element.SerialNumOfIndex(i, j, n);
			k5 = element.SerialNumOfIndex(i + num, j, n);
			k6 = element.SerialNumOfIndex(i - num, j + num, n);
			k7 = element.SerialNumOfIndex(i, j + num, n);
			k8 = element.SerialNumOfIndex(i + num, j + num, n);
			y[element.SerialNumOfIndex(i / 2.0, j / 2.0, n / 2)] = 
				(x[k0] + x[k1] + x[k2] + x[k3] + 8.0 * x[k4] + x[k5] + x[k6] + x[k7] + x[k8]) / 16.0;
		}
	}
	return y;
}


template <typename Element = P1Element>
class Injection : public RestrictionOperator<Element>
{
public:
	Injection(size_t _n);
	~Injection() = default;
	virtual const Vector operator*(const Vector& x);
};

TEMPLATE
Injection<Element>::Injection(size_t _n) : RestrictionOperator<Element>(_n) {}

TEMPLATE
const Vector Injection<Element>::operator*(const Vector& x)
{
	size_t n = RestrictionOperator<Element>::n;
	Element element;
	if (x.rows() != element.GetNodeNum(n))
		throw std::length_error("Injection: sizeof(x) error!");
	Vector y(element.GetNodeNum(n / 2));
	y.setZero(element.GetNodeNum(n / 2));
	double num = sqrt(element.GetNodeNum(n)) + 1.0;
	num = n / num;

	for (double i = 2 * num; i <= n - 2 * num; i += 2 * num)
		for (double j = 2 * num; j <= n - 2 * num; j += 2 * num) {
			y[element.SerialNumOfIndex(i / 2.0, j / 2.0, n / 2)] = x[element.SerialNumOfIndex(i, j, n)];
		}
	
	return y;
}


template <typename Element = P1Element>
class InterpolationOperator
{
protected:
	size_t n;

public:
	InterpolationOperator(size_t _n);
	virtual ~InterpolationOperator() = default;
	virtual const Vector operator*(const Vector& x) = 0;
	size_t GetN() const;
};

TEMPLATE
InterpolationOperator<Element>::InterpolationOperator(size_t _n) : n(_n) {}

TEMPLATE
size_t InterpolationOperator<Element>::GetN() const
{
	return n;
}

template <typename Element = P1Element>
class Linear : public InterpolationOperator<Element>
{
public:
	Linear(size_t _n);
	~Linear() = default;
	virtual const Vector operator*(const Vector& x);
};

TEMPLATE
Linear<Element>::Linear(size_t _n) : InterpolationOperator<Element>(_n) {}

TEMPLATE
const Vector Linear<Element>::operator*(const Vector& x)
{
	Element element;
	size_t n = InterpolationOperator<Element>::n;
	if (x.rows() != element.GetNodeNum(n / 2))
		throw std::length_error("Linear: sizeof(x) error!");
	Vector y(element.GetNodeNum(n));
	y.setZero(element.GetNodeNum(n));
	double value;
	double num = sqrt(element.GetNodeNum(n)) + 1.0;
	num = n / num;
	for (double i = 2 * num; i <= n - 2 * num; i += 2 * num) {
		for (double j = 2 * num; j <= n - 2 * num; j += 2 * num) {
			value = x[element.SerialNumOfIndex(i / 2.0, j / 2.0, n / 2)];
			y[element.SerialNumOfIndex(i, j, n)] += value;
			y[element.SerialNumOfIndex(i, j - num, n)] += 0.5 * value;
			y[element.SerialNumOfIndex(i + num, j, n)] += 0.5 * value;
			y[element.SerialNumOfIndex(i, j + num, n)] += 0.5 * value;
			y[element.SerialNumOfIndex(i - num, j, n)] += 0.5 * value;
			y[element.SerialNumOfIndex(i - num, j - num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i + num, j - num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i + num, j + num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i - num, j + num, n)] += 0.25 * value;
		}
	}

	return y;
}



template <typename Element = P1Element>
class Quadratic : public InterpolationOperator<Element>
{
public:
	Quadratic(size_t _n);
	~Quadratic() = default;
	virtual const Vector operator*(const Vector& x);
};

TEMPLATE
Quadratic<Element>::Quadratic(size_t _n) : InterpolationOperator<Element>(_n) {}

TEMPLATE
const Vector Quadratic<Element>::operator*(const Vector& x)
{
	Element element;
	size_t n = InterpolationOperator<Element>::n;
	if (x.rows() != element.GetNodeNum(n / 2))
		throw std::length_error("Quadratic: sizeof(x) error!");
	Vector y(element.GetNodeNum(n));
	y.setZero(element.GetNodeNum(n));
	
	double value;
	double num = sqrt(element.GetNodeNum(n)) + 1.0;
	num = n / num;
	for (double i = 2 * num; i <= n - 2 * num; i += 2 * num) {
		for (double j = 2 * num; j <= n - 2 * num; j += 2 * num) {
			value = x[element.SerialNumOfIndex(i / 2.0, j / 2.0, n / 2)];
			y[element.SerialNumOfIndex(i, j, n)] += value;
			y[element.SerialNumOfIndex(i - num, j - num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i + num, j - num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i + num, j + num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i - num, j + num, n)] += 0.25 * value;
			y[element.SerialNumOfIndex(i, j - num, n)] += 0.5625 * value;
			y[element.SerialNumOfIndex(i + num, j, n)] += 0.5625 * value;
			y[element.SerialNumOfIndex(i, j + num, n)] += 0.5625 * value;
			y[element.SerialNumOfIndex(i - num, j, n)] += 0.5625 * value;
			if (j - 3 * num >= num)
				y[element.SerialNumOfIndex(i, j - 3 * num, n)] -= 0.0625 * value;
			if (i + 3 * num <= n - num)
				y[element.SerialNumOfIndex(i + 3 * num, j, n)] -= 0.0625 * value;
			if (j + 3  * num <= n - num)
				y[element.SerialNumOfIndex(i, j + 3 * num, n)] -= 0.0625 * value;
			if (i - 3 * num >= num)
				y[element.SerialNumOfIndex(i - 3 * num, j, n)] -= 0.0625 * value;
		}
	}
	return y;
}

#undef TEMPLATE
#endif
